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In cold dark matter cosmological models 1 ' 2 , structures form and grow by merg- 
ing of smaller units 3 . Numerical simulations have shown that such merging is 
incomplete; the inner cores of halos survive and orbit as "subhalos" within their 
hosts 4 ' 5 . Here we report a simulation that resolves such substructure even in 
the very inner regions of the Galactic halo. We find hundreds of very concen- 
trated dark matter clumps surviving near the solar circle, as well as numerous 
cold streams. The simulation reveals the fractal nature of dark matter clustering: 
Isolated halos and subhalos contain the same relative amount of substructure and 
both have cuspy inner density profiles. The inner mass and phase-space densi- 
ties of subhalos match those of recently discovered faint, dark matter-dominated 
dwarf satellite galaxies 6-8 , and the overall amount of substructure can explain 
the anomalous flux ratios seen in strong gravitational lenses 9 ' 10 . Subhalos boost 
gamma-ray production from dark matter annihilation, by factors of 4-15, rela- 
tive to smooth galactic models. Local cosmic ray production is also enhanced, 
typically by a factor 1.4, but by more than a factor of ten in one percent of 
locations lying sufficiently close to a large subhalo. These estimates assume that 
gravitational effects of baryons on dark matter substructure are small. 

The cold dark matter (CDM) model has been remarkably successful at describing the 
large-scale mass distribution of our Universe from the hot Big Bang to the present. However, 
the nature of the dark matter particle is best tested on small scales, where its interaction 
properties manifest themselves by modifying the structure of galaxy halos and their sub- 
structure. CDM theory predicts that the growth of cosmic structures begins early, on Earth 
mass scales 11 ' 12 , and continues from the bottom up until galaxy clusters form that are twenty 
orders of magnitude more massive. Resolving small-scale structures is extremely challenging, 
as the range of lengths, masses, and timescales that need to be simulated is immense. We 
have performed the highest precision calculation - dubbed "Via Lactea II" - of the assembly 
of the Galactic CDM halo. The simulation follows the growth of a Milky Way-size system in a 
ACDM Universe from redshift 104.3 to the present. It provides the most accurate predictions 
on the small scale clustering of dark matter and the first constraints on the local subhalo 
abundance and properties. We used the parallel treecode PKDGRAV2[ 13 ] and sample the 
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galaxy-forming region with 1.1 x 10 9 particles of mass 4, 100 M & . Cosmological parameters 
were taken from WMAP [ 14 ], see the online supplement for more details and a comparison 
to our previous Via Lactea simulation. 

The wealth of substructure that survives the hierarchical assembly process to the 
present epoch is clearly seen in Figure 1: we resolve over 40,000 subhalos within 402 
kpc of the center and they are distributed approximately with equal mass per decade of 
mass over the range 10 6 — 10 9 M Q . They have very high central phase-space densities 
(> 10~ 5 M Q pc~ 3 km~ 3 s 3 ) due to their steep inner density cusps and their relatively small 
internal velocity dispersions. This agrees well with the extremely high phase-space densi- 
ties inferred from stellar motions within ultra faint dwarf galaxies 7 . Our predicted inner 
subhalo densities (0.4 — 2.5 M Q pc~ 3 within 100 pc, 7 — 46 M pc~ 3 within 10 pc) are also 
in excellent agreement with the observations 6,7 . The fact that CDM naturally predicts a 
small-scale dark matter distribution that matches the observations is a real success of the 
model. Particle candidates that introduce a low phase-space limit, such as a sterile neutrino, 
or that have a high collisional cross section such as self interacting dark matter would fail 
these fundamental observational tests. 

The phase-space map (upper inset in Figure 1) also contains coherent elongated fea- 
tures. These are streams which form out of material removed from accreted and disrupted 
subhalos. The few visible streams have quite low densities (about 100 times below the lo- 
cal density) but owing to their low velocity dispersion (about 10 times smaller than that 
of background particles) they just barely manage to stand out in local phase space density 
(these streams have about 10 -9 M pc -3 km -3 s 3 ). These resolved streams together with the 
multitude of expected finer grained phase space structures that we currently can not resolve, 
will lead to unique signatures in direct detection experiments, especially those with direc- 
tional sensitivity. In cases where the disrupted subhalo hosted a luminous satellite galaxy, 
the resulting streams would contain not only dark matter but also stars. This process would 
then produce detectable features in the Milky Way's stellar halo, like those observed in the 
"Field of Streams" 15 . 

Further evidence for halo substructure comes from the anomalous flux ratios in multiply- 
imaged gravitationally lensed quasars 16 ' 17 . Perturbations to the light path from substructure 
can explain this phenomenon if the projected substructure fraction within 10 kpc is about 
one percent 9 ' 10 . Within a projected distance of 10 kpc from the center, 0.50% of the host 
mass belongs to resolved substructure, which could be just enough to explain the observed 
flux anomalies. In earlier simulations this fraction was lower, even the first Via Lactea halo 18 
had only 0.25% indicating that this quantity has not yet converged. 

Via Lactea II predicts a remarkable self-similar pattern of clustering properties: Our 
simulation is the first to use an extremely accurate integration of particle orbits in high 
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density regions 19 , allowing a precise determination of the density profile within the inner 
kpc of the Galactic halo and within the centers of its satellite galaxies. We find that a 
cuspy profile fits the host halos density profile well, while the best fit cored profile lies below 
the simulated inner densities (Figure 2). The inner profiles of subhalos are also consistent 
with cusps over their resolved ranges. They scatter around the moderate cusp index of 
the host halo (7 = 1.24): Some of them are denser in the inner part, and some are less 
dense, exactly like the inner parts of field halos, which have inner slopes of 7 ~ 1.2 ± 0.2 [ 20 ]. 
At large radii subhalo density profiles generally fall off faster than field halo profiles. These 
similarities and differences between subhalo and field halo profiles have a simple explanation: 
Subhalo density profiles were modified by tidal mass loss, which removes material from the 
outside in, but does not change the inner cusp structure 18 ' 21 . Figure 3 shows that the 
dwarf satellites of the Milky Way appear to be scaled versions of the main halo not just 
in their inner mass distribution, but also in term of relative substructure abundances. Via 
Lactea II demonstrates the fractal-like appearance of the dark matter by resolving the second 
generation of surviving sub-substructures from the merging hierarchy. This suggest that at 
infinite resolution one would find a long nested series of halos within halos within halos etc., 
reminiscent of a Russian Matryoshka doll, all the way down the first and smallest earth mass 
haloes that form. 

The multitude of dark substructures increases the dark matter annihilation signal, since 
it is proportional to the square of the local density. For cuspy profiles (Figure 2) with some 
fixed inner slope (7 < 1.5) one gets the following simple scaling relation for the annihilation: 

L OC p 2 s rl OC KtaxAVmax OC V^y/ty , (1) 

see the online supplement for the definition of the concentration cy and its values. Combined 
with the steep subhalo velocity function N(> V max ) oc V~^ x , this implies that subhalos of 
all sizes contribute about equally to the total signal coming from the Galactic dark halo. 
Taking the higher concentrations of smaller systems 22 into account, one finds that small 
subhalos are contributing more than large ones 23 ' 24 . Summing up ^ max /r Vm ax of all resolved 
subhalos in Via Lactea II comes close (97%) to the host halo's ^ax/ r Vmax, i-e. the resolved 
subhalos already contribute as much as their smooth host alone would. In other words the 
"substructure boost factor" is at least 1.97. Extrapolating down to micro-subhalos of size 
Knax = 0.25 ms _1 , taking into account how concentrations depend on subhalo size 22 and 
position (see the online supplement), and assuming a uniform distribution of subhalo inner 
slopes a between 1.0 and 1.5, leads to a total boost of 14.6. Most of it comes from very 
small clumps: halting the same extrapolation at V m£LX = 44ms _1 lowers the boost to 6.6. 
While the contribution from small, dark clumps itself is not affected by baryons, it may 
not dominate the total signal in scenarios in which baryonic collapse greatly increases the 
central dark matter densities in larger halos. However, the net effect of stars, black holes, 
and galaxy formation is unclear, and it may actually lead to a reduction in the central dark 
matter densities. The detailed distribution of cusp indices is still unknown, since only a few 
halos have been simulated with sufficient resolution 20 . For the annihilation boost factors 
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the existence of a few steep cusps near 1.5 would make a big difference, since the signal 
diverges logarithmically towards the center in a 7 = 1.5 cusp. Cutting the assumed uniform 
distribution of inner slopes at 1.4 instead of 1.5 gives a boost of 9.9 instead of 14.6, and 4.3 
instead of 6.6. These factors imply that most of the extra-galactic 7-ray background from 
dark matter annihilation 23 , which will be constrained or even detected by the upcoming 
GLAST mission, should be emitted by subhalos, and not by distinct host halos. 

Besides 7-rays, dark matter annihilation would produce charged particles and anti- 
particles that, due to to magnetic field entanglement, propagate over much smaller dis- 
tances within the Galaxy. Space based experiments (like PAMELA, and in the near fu- 
ture AMS-02) could detect anti-particles produced in dark matter annihilations within 
about one kiloparsec 25 . What fraction of this local annihilation would happen in nearby 
subhalos? To constrain this local boost factor we use the same assumptions as above 
(7 = [1 — 1.5],y max > 0.25 ms -1 ), but now we only include subhalos within one kilopar- 
sec of the solar system (see the online supplement for the local subhalo abundance). The 
resulting signal is 40% of the smooth halo signal, giving a boost factor of 1.4, which we esti- 
mated using the spherically averaged density at 8 kpc (p = 0.40 GeV c~ 2 cm -3 ). Explaining 
the positron excess from the HEAT experiment 26 with local dark matter annihilation requires 
enhancements from about 3 to 100 25 . When a relatively large subhalo happens to lie within 
1 kpc, one can get the higher boost factors without violating the local subhalo constraints 
from our simulation. Such cases are possible, but not likely: Only 5.2 percent of all random 
realizations have a boost factor of 3 or larger (caused by a V max > 3.4 km s" 1 clump within 1 
kpc). In only 1.0 percent of the cases the boost factor reaches 10 or higher due to a nearby, 
large V max > 5.6 km s -1 subhalo. 
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Figure 1: Projected dark matter density-square map of "Via Lactea II". An 800 kpc 

cube is shown. The insets focus on an inner 40 kpc cube, in local density (bottom), and in local phase 
space density calculated with EnBiD[ 27 ] (top). The Via Lactea II simulation has a mass resolution 
of 4,100 M© and a force resolution of 40 pc. It used over a million processor hours on the "Jaguar" 
Cray XT3 supercomputer at the Oak Ridge National Laboratory. A new method was employed 
to assign physical, adaptive time-steps 19 equal to 1/16 of the local dynamical timescale (but not 
shorter than 268,000 yr), which allows to resolve very high density regions. Initial conditions were 
generated with a modified, parallel version of GRAFIC2[ 28 ] . The high resolution region is embedded 
within a large periodic box (40 comoving Mpc) to account for the large scale tidal forces. The mass 
within 402 kpc (the radius enclosing 200 times the mean matter density) is 1.9 x 10 12 Mq. 




Figure 2: Density profiles of main halo and subhalos. Main panel: Profile of the 
Milky Way halo (thick line) and of eight large subhalos (thin lines). The lower panel gives the 
relative differences between the simulated main halo profile and a fitting formula with a core 29 
p(r) = /9 s exp{— 2/a [(r/r s ) a — 1], with best fit parameters: a = 0.170, r s = 21.5 kpc, p s = 1.73 x 
10~ 3 Mq pc~ 3 (red curve) and one with a cusp 20 p(r) = p s (r /r s )~' ) '(r jr s + 1) _3+7 with a best fit 
inner slope of 7 = 1.24, r s = 28.1 kpc, p s = 3.50 x 10~ 3 Mopc~ 3 (blue curve). The vertical dotted 
line indicates the estimated convergence radius of 380 pc: simulated local densities are only lower 
limits inside of 380 pc and they should be correct to within 10% outside this region. The cuspy 
profile is a good fit to the inner halo, while the cored profile has a too shallow slope in the inner 
few kpc, causing it to overestimate densities around 4 kpc and to underestimate them at all radii 
smaller than 1 kpc. The simulated densities are higher than the best fit cored profile even at 80 pc, 
where they are certainly underestimated due to numerical limitations. We find the same behavior 
in the inner few kpc in all six snapshots we have analyzed so far between z=3 an z=0. The large 
residuals in the outer halos on the other hand are transient features, they are different in every 
snapshot. Inset: Rescaled host (thick line) and subhalo (thin lines) density profiles multiplied by 
radius square to reduce the vertical range of the figure. 




Figure 3: Subhalo and sub-subhalo abundances. Number of subhalos above V max 
within r2oo — 402 kpc (thick solid lines) and within 100 and 50 kpc of the galactic center 
(thin solid lines). V m3X is the peak height of the subhalo circular velocity v C i rc = JGM(< r)/r 
and serves as a simple proxy for the mass of a subhalo. The dotted line is iV(> V max ) = 
0.036 (Vmax/Vmax^ost) -3 , where Vma^host = 201 km s" 1 (at r V max,host = 60 kpc). It fits the 
subhalo abundance above V^ax ~ 3.5 km s -1 . The number of smaller subhalos is artificially 
reduced by numerical limitations. Inside r2oo this halo has 1.7 times more substructure 
than the first Via Lactea halo 18 , a factor well within the halo-to-halo scatter 30 . Inside 50 
kpc the difference grows to 2.6, probably due to the improved mass and time resolution 
of Via Lactea II, which allows to resolve inner substructure better. The inset shows the 
sub-subhalo abundance within r 10 oo (enclosing 1000 times the mean matter density) of the 
centers of eight (same ones as in Fig. 2) large subhalos (thin solid lines), riooo is well inside 
of the tidal radius for these systems. The thick solid line shows the subhalo abundace of 
the host halo inside of its riooo = 213 kpc. The (sub-)subhalo Vm ax values are given in 
units of Viooo = \JGM{< riooo) /^iooo of the corresponding host (sub-)halo. Lines stop at 
Vmax = 2kms~ 1 . The mean sub-substructure abundance is consistent with the scaled down 
version of main halo, and both the mean abundance and the scatter agree with results in 30 
for distinct field halos. 
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In this Supplementary Information, we give more details about the Via Lactea II (VL- 
II) simulation presented in our Letter. The first part provides additional results about radial 
trends in subhalo abundance and properties. In the second part we include some comparisons 
with the Via Lactea (VL-I) run, previously the largest simulation of the Galactic dark matter 
halo. 

Supplementary Figure 1 illustrates how tidal interactions with the gravitational poten- 
tial of the host shape the radial distribution of subhalos and their internal structure. Tides 
remove mass from the outer parts of subhalos. Clumps closer to the Galactic center lose more 
mass as they experience stronger tidal forces and more pericenter passages 18 . Mass loss also 
reduces V max and r Vm ax, which leads to an increase in p{< r Vm ax)- The radial distribution 
of subhalos with Knax > 3 km s" 1 is more extended than the dark matter distribution in the 
Galactic halo, a feature that does not depend on subhalo size, i.e. different V mSuX selection 
thresholds lead to the same radial distribution. VL-II resolves subhalos as close as 8 kpc from 
the Galactic center, but it is possible that it underestimates the true substructure abundance 
inside 20 kpc due to numerical limitations. Subhalo concentrations are defined as the mean 
density within ry max , the radius of peak circular velocity, a quantity that is well determined 
both for isolated halos and subhalos and does not dependent on assumptions about their 
density profiles 18 : 

Cy = p(< r Vmax )/Pcrit,0 = (KnaxAVmax) 2 ^ / (3Gp C rit,o) , (2) 

where p C rit,o — 1-48 x 10~ 7 M pc~ 3 . The median subhalo concentration increases strongly 
towards the Galactic center, both because of tidal mass losses and to a lesser extent because 
of the earlier formation times of inner substructure 18 . 

Supplementary Table 1 summarises the numerical parameters and host halo properties 
of the VL-I and VL-II simulations. The main difference between these runs is the improved 
time stepping in VL-II, which at each time and for each particle is based on the local 
dynamical time ^ 'l / G p enc [ 19 ] , where p enc is the mean density enclosed in a sphere extending 
out to the particle's position and centered on the dynamically dominant structure. In VL-I 
we adopted instead the standard ad-hoc time-step criterion based on the acceleration and 
the gravitational softening of the particle (oc y / e/jaj). This criterion fails in high density 
regions and artificially flattens the inner cuspy halo density profiles 31 . This limitation set 
the convergence radius of VL-I (where the host true density profile is reproduced to within 
10 percent) to r conv = 1.3 kpc. The new time stepping used in VL-II allows us to properly 
resolve the host density profile on significantly smaller scales: the finite mass resolution 
determines a convergence radius 20 of about 0.38 kpc for VL-II. The VL-I subhalo density 
profiles are also affected by this limitation: the enclosed densities within 300 pc of large 
subhalos are about twice as high in VL-II. Further out, at 600 pc, the enclosed subhalo 
densities are very similar in VL-I and VL-II. 
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To check for numerical convergence and test the dependence of our results on numerical 
and cosmological parameters we ran a series of lower resolution versions of VL-II. The mass 
resolution in this "VL-IIm" series is 64 times coarser and the force softening length 4 times 
larger than in VL-II. Supplementary Figure 2 shows that the lower resolution version of the 
same halo has a very similar subhalo velocity function above about 0.05 Vmax,host — 10 km s _1 . 
Rescaling to 64 times less massive systems suggests that VL-II should have converged down 
to about 2.5 kms^ 1 , which indeed is close to the scale where the VL-II velocity function 
starts to fall below the power law fit (Figure 3). The earlier starting redshift of VL-II does 
not seem to affect the z = substructure abundance significantly. We find only a weak 
dependence on cosmological parameters: VL-II used the best fit ACDM parameters from 
the WMAP 3 year data release 14 : tt m = 0.238, fi A = 0.762, h=0.73, n s =0.951, and cr 8 =0.74. 
For comparison we have ran a simulation with WMAP 1 year parameters 32 : Q m = 0.27, 
Q\ = 0.73, h=0.72, n s =1.0, and er 8 =0.9. The higher cr s and steeper spectral index n s in 
WMAP1 lead to more small scale power. The effect on the z = substructure abundance 
is rather small: our WMAP1 run has about 20 to 30 percent more subhalos relative to the 
WMAP3 runs, in agreement with semi- analytical predictions (see Figure 11 in 33 ). 

31. Diemand, J., Zemp, M., Stadel, J., Moore, B. & Carollo, C. M. Cusps in cold dark 
matter haloes. Mon. Not. R. Astron. Soc. 364, 665-673 (2005). 

32. Spergel, D. N. et al. First- Year Wilkinson Microwave Anisotropy Probe (WMAP) 
Observations: Determination of Cosmological Parameters. Astrophys. J. Supp. 148, 175- 
194 (2003). 

33. Zentner, A. R. & Bullock, J. S. Halo Substructure and the Power Spectrum. 
Astrophys. J. 598, 49-72 (2003). 
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Supplementary Figure 4: Abundance and concentrations of subhalos vs. distance 
from the galactic center. Top: The number density profile of subhalos (circles) is more 
extended than the dark matter density profile p(r) (thick line). Their ratio turns out to be 
roughly proportional to the enclosed mass M(< r), i.e. pM(< r) (thin line) matches the 
subhalo number density quite well. Only subhalos larger than V max = 3 km s" 1 are included 
here. Bottom: Subhalo concentrations (median and 68% range are shown) increase towards 
the center, where the stronger tidal force remove more of the outer, low density parts from 
the subhalos. To make sure their cy are resolved, only subhalos larger than V mauX = 5kms _1 
are used. The error bars indicate the statistical uncertainties in both panels. 
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Supplementary Figure 5: Subhalo abundance at different numerical resolutions, 
starting redshifts and cosmologies. Number of subhalos above Knax/Knax^ost within 
r 20 o for the VL-II simulation and three lower resolution versions of the same halo. 
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Supplementary Table 1: Simulation parameters and halo properties. Time step 
criterion AT, spline force softening length e, initial redshift z iy convergence radius of the 
host halo density profile r conv , mass M hires of high resolution dark matter particles and host 
halo r 20 o, M 2 oo, V mSuX and rVmax for the VL-I and VL-II simulations. At each time individual 
particle time steps are chosen by dividing the base time step of 13.7 Gyr / 400 by two until it 
is smaller than AT, where \a\ is the norm of the acceleration vector and p cnc is the enclosed 
density within the dynamically dominant structure. Force softening lengths e are constant 
in physical units back to z = 9 and constant in comoving units before. 
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